# Perform Gausian Kernel density approximation
# Ref to: Point, 2014, pp. 83-84
library("spatstat")

data <- read.csv("input/data.csv")

# calculate sigma optimums in 2 different ways
sigma <- (sd(data$zinc) + sd(data$copper))/2
# way #1
sigma.optim1 <- 1.06 * sigma * nrow(data)^-0.2
# way 2
iqr <- (IQR(data$zinc) + IQR(data$copper))/2
sigma.optim2 <- 0.9 * min(sigma,iqr)*nrow(data)^-0.2


# create object with 2 Z-fields: zinc & copper
object <- ppp(x = data$x, y = data$y, marks = data.frame(zinc = data$zinc, copper = data$copper), window = owin(c(min(data$x),max(data$x)), c(min(data$y), max(data$y))), unitname = "meters")
plot.ppp(object, use.marks = T)

fit <- density.ppp(object, sigma.optim2, edge = T, kernel = "gaussian")
image(fit, which.marks = "zinc")
plot(object, which.marks = "zinc", add = T)
